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ABSTRACT 

Gravitating bodies significantly alter the flow pattern (density and velocity) of the gas 
that attempts to stream past. Still, small protoplancts in the Mars-super-Earth range 
can only bind limited amounts of nebular gas; until the so-called critical core mass 
has been reached (~1-10 Earth masses) this gas is in near hydrostatic equilibrium 
with the nebula. Here we aim for a general description of the flow pattern surrounding 
these low-mass, embedded planets. Using various simplifying assumptions (subsonic, 
2D, inviscid flow, etc), we reduce the problem to a partial differential equation that we 
solve numerically as well as approximate analytically. It is found that the boundary 
between the atmosphere and the nebula gas strongly depends on the value of the 
disc headwind (deviation from Keplerian rotation). With increasing headwind the 
atmosphere decreases in size and also becomes more asymmetrical. Using the derived 
flow pattern for the gas, trajectories of small solid particles, which experience both 
gas drag and gravitational forces, are integrated numerically. Accretion rates for small 
particles (dust) are found to be low, as they closely follow the streamlines, which curl 
away from the planet. However, pebble-size particles achieve large accretion rates, in 
agreement with previous numerical and analytical works. 

Key words: planets and satellites: formation - planetary systems: protoplanetary 
discs - hydrodynamics - methods: analytical 



1 INTRODUCTION 

Accretion of gas by gravitating objects is a common theme 
in astrophysi cs. The pioneering work i n this field goes back 
the w orks of lHovle fc Lvttletonl (|l939h and lBondi fc Hovld 
1 19441) . who analytically derived the accretion rates for 
point-bodies (together known as Bondi-Hoyle-Lyttleton or 
BHL-accretion) . Follow-up studies have refined these esti- 
mates, by considering the effects of the Mach number of the 
incident flow and the equation of state for the gas, mostly for 
super sonic flows (e.g., lHuntlll979l ; IShima et al.|[l985l : lRuffertl 
Il995l ). In astrophysics, BHL-accretion plays a role in Binary 
systems, star- and ga laxy formation, and the formation of 
gas-rich planets (e.g . , iKIev et al.1 Il995l : iBonnell etail l200ll : 
iNelson fc Benzll2003l ). 

According to the core accretion model of planet forma- 
tion, pla nets start out small and rocky a s protoplanetary 



embr yos (lMizundll980l; iPollack et al.lll996l ; IWeidenschillinel 
ll997l :' lHubickvi et " al.ll2005l )7 Once these embryos reach a size 



where their Bondi radius 



Rb 



G N M P 



(1) 



the gas from the nebula, forming a putative atmosphere. In 
equation |T|) Gn is Newton's constant, M v the mass of the 
perturber, and Coo the sound speed of the unperturbed sys- 
tem. As accreting planets are hotter than the background 
nebula - either from infalling solids ( planetesimals) or fro m 
the residual heat of their formation (llkoma fc Horill2012l ) - 
the resulting pressure gradient prevents contraction of the 
atmosphere. Therefore, the amount of gas these (small) bod- 
ies can acquire is initially limited: the protoplanet's atmo- 
sphere is in pressure-equilibrium with the surrounding neb- 
ular gas. Small, hot protoplanets do not acquire gas akin 
to the BHL-regime: they must cool first. 

However, at some critical core mass M cr it the steady 
state picture is no longer appropriate. Many detailed, ID 
models have been developed to describe the quasi-static den- 
sity structu re of the protop l anet atmosphere and the cross- 
over point (IStevensonl l 19821 : IWuchterl|[l993l ; llnaba fc Ikomal 
120031 : lRafikovll200(i r This crossover mass is often quoted 
as ~10 M 9 (Earth masses) but its precise value depends 
on detailed atmosphere models, and thereby on the equa- 
tion of state, accretion luminosity, and opacities of gas and 



starts to exceed their physical radius R, they start to bind 
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1 Ind eed, in a numerical experiment iTerquem fc H cincmami 
ll201lh recently showed that, starting from a much condensed con- 
figuration, a protoplanet's atmosphere will expands and reconnect 
with the disc's gas to restore the near-hydrostatic state. 
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grains; it decreases to ~1 for grain- free atmospheres 



|Hori fc Ikomal I2010T ). Obviously, understanding the pro- 
cesses that determine Merit is of key importance for planet 
formation; it may determine the transition between (super)- 
Earth/Neptune-like planets (bodies that have acquired lim- 
ited amounts of gas) and gas giants (bodies that accrete gas 
in the BHL-regime). 

In this work we will focus on the small planet regime, 
assuming a quasi-steady state. In a disc, the unperturbed gas 
flow, as seen from a frame rotating with the planet, consist 
of two components: a shear (due to the [Keplerian] rotation 
of the disc) and a systematic offset, which we refer to as the 
headwind. The latter arises from the slightly subkeplerian 
rotation of the gas, which is partially support by pressure. 
For a planet o n a circular orbit, the value of the headw ind 
is rjv k with 77 (jAdachi et al.lll97d ; IWeidenschillindll977t ): 



'/ : 



dP 



2paQ? da 



1 c 2 3 dlnP 

2 vj( d In a 



(->) 



where a, p, P, fi, and vk are the disc radius (semi-major 
axis), the corresponding gas density, pressure, orbital fre- 
quency, and orbital velocity. In this work the headwind is 
also defined in terms of a Mach number, «h w = -MhwCoo, i.e., 



Mi 



r)v k /c 



n l/2 



Only positive values of «h w (a head- 



wind; the gas rotates lower than Keplerian) are consider, but 
we remark that in special locations ('pressure bumps') «h w 
may reverse sign. In addition, the flow, as seen in the frame 
of the planet, is subject to (noninertial) Coriolis forces and 
to the tidal force from the distant star. The flow pattern is 
therefore quite rich. 

Understanding the flow pattern close to the planet is 
important for several reasons. For example, the flow pat- 
tern in the co-orbital region of the planet, where stream- 
lines make a 'U-turn' (the horses hoe) is a c ritical in- 
gredient for the co-orbital torq ue (|Wardl Il99ll ; see also 
iPaardekooper fc Papaloizoull2009l ). Also, the location of the 
boundary between the nebular gas and the protoplanet's at- 
mosphere is a parameter that affects the therm al evolution of 
(growing) protoplanets (|Lissauer et al.l [20091 ). The amount 
of gas that can be bound during the protoplanetary disc 
phase affects t he long-term evolution of super-Earths and 
mini-Neptunes lllkoma fc Gendal 120061 ; llkoma fc Horl |2012| ; 
lLopez et alj|2012ft . 

The gravitationally perturbe d flow pattern also af- 
fects the b ehavior of pa rt icles. IWeidenschilling fc Davia 
(|l985l ') and IPaardekooper] (|2007f > investigated the accre- 
tion potential of (small) particles. The secular particle- 
planet interaction has been studied more generally by 
iMuto fc Inutsukal (|2009l 'l. Of special importance in the con- 
text of this work is the new 'pebble-accretion' mechanism, 
where protoplanets quickly accrete small particles, as found 
by both numerical simulations (|Lambrechts fc Johanser] 
l2012i ; lMorbidelli fc Nesvornvl2012i) as well as analytical est i- 
1ft I 



mates (|Ormel fc Klahjl201ot iPerets fc Murrav-Clavl l201lT ) . 
The latter studies however did not account for the modifica- 
tion of the gas flow by the gravity of the protoplanet - the 
topic of this paper. 

The goal of this paper is to obtain a quantitative de- 
scription of the density structure and the flow (gas velocity) 
in the vicinity of a gravitating body - in particularly, a pro- 
toplanet in a Keplerian-rotating disc - and to assess the 
role of key parameters like: the headwind, shear, and planet 



mass. The low-mass planet regime under consideration is 
characterized by the following scale hierarchy: 



R < R b < Rh < H, 



(•3) 



where H is the scaleheight of the disc, R the radius of per- 
turber, and Rh the Hill radius 



Rh = a 1 



3M* 



1/3 



(4) 



(with the stellar mass) - the scale on which the so- 
lar gravity rivals that of the perturber. At the lower range, 



R ~ Rb, bodies are about ~ 10 km in size, or 10 



-1(T 



Mgj in mass. With increasing mass both Rb and Rh increase, 
but the ratio Rh/Rb decreases. At the high-mass end of in- 
equality Q Rb — Rh — H. This corresponds to masses of 
~10 M®, somewhat larger for the outer disc. 

One key approximation that is employed is to neglect 
contributions originating from the Lindblad region, i.e., the 
distance x ~ H, where shock waves are excited. Resolv- 
ing the shock at these distances is important for calculating 
the torque that is exerted on the planets, which determines 
the orbital evolution (migration) of the planet. However, as 
an approxi mation for the flow close to the pla nets, i.e., on 
scales <C-ff. IPaardekooper fc Papaloizoul l|2009l 'l argued that 
it is justifiable to omit the contributions arising from the 
Lindblad torque region. That is, they showed that this pro- 
cedure, referred to as 'torque cut-off', is permissible - up 
to levels of ~10% accuracy - provided that the softening 
radius, which is customary applied to smoothen the gravi- 
tational potential in numerical simulations, is chosen small 
enough. In this study no softening is present; instead the 
numerical setup is characterized by a surface (inner bound- 
ary condition). Avoiding the softening is advantageous, since 
it is a parameter that, unless carefully ch osen, could affect 
the outcome of the numerical experiment (|Dong et al.ll201ll : 
iMiiller et aHl2012h . 

These considerations imply that the flow past the pro- 
toplanet is subsonic. A further assumption is that the flow 
is two dimensional (2D), which allows us to formulate the 
flow in terms of a single scalar quantity, the stream function 
^t. Employing the 2D assumptio n significantly reduces the 
comp lexity of the problem (e.g., iKorvcanskv fc Papaloizoul 
Il996f ). Extending the stream function formulation to 3D 
configurations is str aightforward for axisymmetric flows 
(|Lee fc Stahlerll201lf) . but much more difficult - often un- 
practical - for truly 3D flows. Similar simplifying assump- 
tions are that the flow is inviscid and steady (non-turbulent). 
These idealizations allow us to conduct a thorough parame- 
ter study for the flow pattern in the vicinity of the perturb- 
ing body, making a systematic investigating of the sensitiv- 
ity of the flow pattern to the various parameters involved 
(protoplanet mass, nebular headwind, equation of state, nu- 
merical parameters) possible. 

The structure of the paper is as follows. In Section[3]we 
formulate the problem and the underlying (inviscid) fluid 
equations. In Section[3]a first analytic model for the flow pat- 
tern is presented by a linear perturbation analysis. Section 5] 
highlights the key results obtained from a numerical pa- 
rameter study. Section [5] a more complete analytic model 
is given, which describes the critical atmosphere region, un- 
der the approximation that the density is a radial function. 
In Section [6] the analytic model for the gas flow is used to 
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numerically integrate trajectories of small particles in order 
to assess their accretion behavior. Further implications are 
discussed in Section [7] Section [8] presents the conclusions. 



2 STREAM FUNCTION FORMULATION 

In this section the equations of continuity and force-balance 
(Euler) are rewritten in terms of Bernoulli's equation and 
a diffusion equation for the stream function ^ - the quan- 
tity that together with the surface density E characterizes 
the flow. The numerical solution to these equations is pre- 
sented in Section[4]and an analytical approximation is given 
in Section [5] A key feature of the approach is to connect the 
local solution near the perturber to the unperturbed solu- 
tions at large radii - the far field. The far field flow pattern 
is assumed to be known and specified by the following con- 
stant quantities: surface density (£«>), sound speed (coo), 
vorticity (a)c»), and headwind (w h „). 



2.1 The stream function and far field solution 

The 2D shearing-sheet approximation is adopted, with the 
centre of the coordinate frame rotating at an angular fre- 
quency equal to the local orbital frequency $lo- In the ab- 
sence of a perturbing body the flow is, without loss of gen- 
erality, assumed to be directed in the negative y-direction: 



= (— «hw — q&ox)e y = (—Vhw + w oa x)t 



(5) 



where q is the dependence of fl on disc radius, f2 oc a~ q at 
a = ao. All subscripts involving 'oo' refer to the far field 
(unperturbed) solution. Instead of q the linear shear can be 
quantified by the vorticity in the far field: 



: Move, = V X V c 



(6) 



In a Kcplerian disc, q = 3/2 and io~o = — 3£7o/2. 

The 2D stream function, — ^e z , is defined such that 
it satisfies the continuity equation in steady state, V-Eu = 0: 

SiiEVx*e z . (7) 

In the far field, E = Eoo is constant and >t becomes 



*oo = ( v bw x — -uj x x ) E 



(8) 



Dividing the continuity equation Q by E and taking the 
curl we obtain the vorticity of the flow: 

/ V X *\ 1 
W = e 2 .Vx» = e 2 - Vx— U-V-(-V*), (9) 

where we used the (general) vector identity e z • V X A — 
— V ■ (e z X A) and the 2D-specific 



V* = e- X (V X *e z ). 



(10) 



Furthermore, the vortensity (2f2o + u))/T,, is conserved along 
streamlineslfl Since E and u> are constant in the far field, the 



2 See e.g.. lLandau fc Lifshitd i 19591 ) for a proof. Morover, con- 
servation of vortensity implies that the flow is barotropic (equa- 
tion [IS). 



vortensity is constant everywhere and uj can be expressed in 
terms of E: 



u) = (2fio + ^oo)^ 2f2 ; 

and inserted in equation ([9]) to obtain 



-w— + 2fi , 



(11) 



(12) 



where ui = 2Qq + uioo (= +S7o/2 in a Keplerian disc). Thus, 
^ obeys a diffusion equation in which both the diffusion co- 
efficient (1/E) and the source term (the RHS of equation ll2|l 
are functions of E. If the background flow is shear- free (head- 
wind only) the source term is zero. If in addition the flow is 
incompressible Laplace equation is obtained: V \P = 0. For 
the general case, however, both ^ and E are unknown. 

2.2 Bernoulli's equation 

To solve the degeneracy between E and SP, equation (|12[) is 
supplemented by the Euler equation (force balance). In the 
uniformly rotating reference system that is considered here, 
the Euler equation contains fiducial forces: 



dv 

— +v-Vv 

at 



-2n x«-n ( ) x (n xr)--vp-v$ s , (13) 

2-i 



where $ 9 is the combined gravitational potential (of the 
planet and the star) acting on the fluid element. The pres- 
sure gradient term VP/E is split into a global contribu- 
tion, which gives rise to the headwind, and a local per- 
turbation caused by the planet. The global term equals 
(c(P 00 /dr)/E 00 = — 2?7«fcf2o = — 2vh w fio (see equation[2]) and 
is added to equation (|13[) such that VP/E denotes the per- 
turbation. In addition, it can be shown that the centrifugal 
force term (— flo X (flo X r) — — fio'" 2 ) and the solar force 
term (V"3? SU n) combine into 2qClo X e x . Since by assumption 
the flow is steady (d/dt — 0), equation 1)131) reads 

1 i 

v ■ Vd = -2U xv- —VP + {2qQ, 2 x + 2v hv {l )e x - V$p, 

(14) 

where <E>p is solely due to the gravitating body at the centre 
of the reference frame and VP corresponds to the pressure 
perturbations induced by it. In the far field the LHS as well 
as the pressure and potential terms on the RHS of equa- 
tion (|14|) vanish; and it can be verified that equation ([5j is 
a solution to equation (| 14[) . 

Using the identity \Vv 2 = v ■ V« + v X V X v equa- 
tion (|14[) transforms into: 

1 9 1 

2ft Xv-v XV Xv =- -Vw - -VP (15) 

- V ($p - qQ, x 2 - 2v hv ,il x) , 
or, in terms of the vorticity lj = V X v. 
v X (w + 2Q e z ) = VB (16) 
with B the Bernoulli function 

B = ^v 2 + W + $p + w^Qox 2 - 2v hw n x. (17) 



In equation (|16[) it was further assumed that the flow is 
barotropic: 



(18) 
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Dimensionless mass m = Rb/H = q P /h 3 



Figure 1. Left j/-axis: the dimensionless Hill radius rh as func- 
tion of mass parameter m = R^/H (solid line). Right y-axis: the 
perturber mass M p corresponding to m for two values of the disk 
aspect ratio h = H/ao (dashed lines). For the conversion of q p 
to M p the star was assumed to be of solar mass. Vertical dotted 
lines show the mass where R = Rj,, the point where the embryo 
starts to bind the nebular gas, at distances of 1 AU and 10 AU. 
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Figure 2. (gray curves) The enthalpy function W (equation 128 jl 
as function of density for several values of the transition den- 
sity: = 10, 10 3 and oo (isothermal), (black curve) the LHS of 
equation (25} for a value for ± | V*| 2 = 10 4 . The minimum density 
/min for which the Bernoulli equation will obtain real solutions is 
indicated by an open circle. 



with W the enthalpy. In the far field W is assumed to be 
zero, VK(Eoo) = 0, and the Bernouilli function reduces to 



-Boo = -«» + WoattoX 2 ~ 2V\ m Q.0X 



(19) 



1 2 



1 2 



where equation ((5JI has been inserted for v x . 

Since the LHS of equation (|16l) is orthogonal to v, it 
vanishes upon integration along a streamline. This implies 
that B is conserved along streamlines (Bernoulli equation). 
To find the relation between B and ^, equation ([7]) is in- 
serted into equation (|16[1 , The LHS of equation (|16|l is then 
written in terms of 



V X *e z S _ w _ T _„ 
X — »e z = -—V* = VB, 

Zj 2-iqo oo 



(20) 



where equation (| 11 f) and identity (|lUp are used. As ui and 
Eoo are constant, equation (|20|l simply implies that B is a 
linear function of 



& T 1 2 



(21) 



where the integration constant has been defined such that 
S(*oo) = -Boo (equation IT9j). 

Together, equations (|17|) , (I21|l provide a second relation 
between *]/ and E: 

1JWJ^_«| 5L+ w y + w+ $ p+WooQox *_ 2vhv/ n oX = 

A 2-1 _■ Z_J rx3 

(22) 



unit, which implies that the unit of velocity is Coo- The di- 
mensionless disk headwind TVfhw = ^hw/c s is referred to as 
the Mach number of the headwind. Because Ri, scales with 
the perturber mass, the dimensionless Bondi unit is also a 
dimensionless mass: 



H 



GM P 



h 3 



(23) 



where q p = Mp/M* and h = H/ao the aspect ratio of the 
disk. We will use m as the parameter for the mass of the per- 
turbing body. The dimensionless Hill radius can be written 



Th 



Rh = I 
H ~ \3h 3 J 



\ 1/3 



= (m/3) 1/3 . 



(24) 



Figure [T] shows the relation among these quantities as func- 
tion of the mass parameter m. The regime m < 1 applies 
in this paper (low mass limit). Small bodies start to bind 
an atmosphere when their physical radius starts to exceed 
the Bondi radius, R ~ 7?;,. This correspond to an m-value of 
m ~ \J M+l p a a\ where p s is the internal density of the body. 
With increasing mass the Bondi and Hill radius start to con- 
verge on each other. The point m ~ 1 signifies the transition 
to the high mass regime (for which 7?;, > Rh > H) which is 
not considered here. 

In nondimensional units the diffusion equation (|12|l and 
Bernoulli equation (I22|l read: 



V- ^V* = -wE + 2 



(25) 



2.3 Units and nondimensionalisation 

For our problem it is natural to measure frequencies ((jJ,lo 
and £1q) in terms of £lo and densities in terms of Eoo- A 
natural length unit is not as obvious, since both the Bondi 
radius, the Hill radius, and the scaleheight of the disk H 
could qualify. Here we adopt the scaleheight as the length 



- l -^y-+W(X) = -Cj*+- + J -^-w oc x 2 +2M h „x. (26) 

These are two equations for the two unknown quantities, E 
and The unperturbed solution (E = Eoo = 1; * = *oo) 
is obtained when m = 0, as it was specified that W(E = 
1) = 0. 
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Figure 3. The perturbed stream function ip y (black solid curve) 
and the radial velocity —v x as function of y in the linear approx- 
imation of the linear model for a mass parameter of m = 10 — 2 . 
The perturbed surface density which follows from equation i31l 
is given by the red thick line. The red thin line gives cr in the 
hydrostatic limit: a = m/r. 



2.4 Equation of state (EOS) 

These equations are supplemented by an equation of state 
(EOS). We adopt: 



P = 



Sx 1 



(E < E T ) 
(E > E T ) 



(27) 



that is, the EOS changes from isothermal to adiabatic 
at a transition density Et- This prescription approxi- 
mat es the structure of embedded protoplanets i n disks 
fsee llnaba fc Ikomall20o3 : lQrmel fc Kobavashj|2012ft . Thus, 
Et is a proxy of the thermodynamic state of the atmo- 
sphere, signifying the point where cooling becomes ineffec- 
tive. Larger (grain) opacities and more luminous protoplan- 
ets will have a smaller isothermal layer, meaning a lower 
Et • 

The enthalpy, through the definition in equation (|18p . 
becomes (in units where = 1): 



W = 



logE 

7 

7-1 



_E_ 

e7 



7-1 



log Et 



(E ^ E T ) 
(E > E T ) 



(28) 

This equation is plotted as function of E in Fig.[5]for several 
transition densities. In equation (|28|) the integration con- 
stants have been chosen such that W is continuous at Et 
and W(0) = 0. The adiabatic index is fixed at 7 = 1.4. 



3 A LINEAR, ANALYTIC SOLUTION TO THE 
PERTURBED FLOW PATTERN 

Assuming that the perturbations are small, equations (|25p . 
(|26[) can be linearized by setting E = 1 + a, — 'foe + ip 
where a and ip are assumed to be small with respect to 
the unperturbed quantities. Inserting these expressions into 
equations (|25l) . (|26[) and keeping only the first order pertur- 
bations (that is, ignoring terms of order higher in ip and cr), 



one obtains (see Appendix [Al for details): 
X7 2 lp + (d x a)(Mhw — w„oi) + er^oo + £>) 



(29) 



cr[l-(o; 00 x-X hw ) 2 ] + (Lj 00 a;-yV(hw)(9 :l; V) = -Qip+ — (30) 

This set of equations, although simpler than the original, 
still represents a complex partial differential equation (PDE) 
which generally can be solved only numerically. To neverthe- 
less pursue with an analytical model, two approximations 
are made: 

(i) The perturbed flow depends only on y, ip = ip(y); 

(ii) The limit 1 < 1 and x <ZCy are considered. 

The first approximation implies that the perturbed flow only 
describes the x component of the velocity. This is reasonable 
since the y-component of the flow is in any case dominated 
by the unperturbed solution (ojoaxey), whereas v x is entirely 
determined by the perturbed component. Thus, d x ip = and 
Vx = (d y ip)/Y< « dyip in the perturbed limit. The second as- 
sumption implies that the focus lies on the co-orbital region, 
x ~ 0. This is again reasonable as is small here. Since 
terms including x axe now small they can be neglected and 
Equation (|3UI) becomes simply 



_ , m 

cr w — uiip H , 

r 



(31) 



where we have also neglected the -M hw term. 

With the same argumentation, the middle term in equa- 
tion (I29|l can be neglected; inserting equation (|31[) into equa- 
tion (|29p and substituting y for r (in fact setting x = 0) 
results in an ordinary differential equation (ODE) for ip(y): 



d^ip — 2<Zi(l + tJoo)^ 



2(1 



(32) 



The solution to this ODE is denoted ip y o. Closed- form so- 
lution exists, but depend strongly on the value of . In a 
Keplerian disk, cjoo = —3/2 it reads: 



IpyO 



[^\y\- 2y Si(y)] cos y + 2y Ci(y) sin y 
V2\y\ 



(33) 



where y — y/y/2 and Ci(u), Si(u) the sine and cosine integral 
functions, defined as 



Ci(u) 



cos(t) 

t 



dt; 



Si(«) 



u shift) 



dt. 



(34) 



The general solution to equation (|32l) contains two integra- 
tion constants. These are of sinusoidal nature, ip ~ smy, 
related to the solutions of the homogeneous equation. These 
oscillatory terms are discarded however, since such oscilla- 
tory terms do not decay for y — > 00. Thus, equation (|33[) 
only contains the particular solution to equation (|32[) . Note 
that equation (|33|l is linear in m. 

Equation (|33|1 . plotted in Fig. [3] peaks at y — for 
which ip v o — iphs = irm/\/2. With the caveat that a is no 
longer small here and that as a consequence the solution will 
break down when y — > 0, this value of ip corresponds, in the 
shear-only case (Alhw = 0), to the widest streamline that 
undergoes so called horseshoe motion - i.e., where the gas 
elements upon their approach to the central object make a 
U-turn. This critical streamline therefore defines the width 
of the horseshoe region (xhs), which, as already alluded to 
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-1.0 -0.5 0.0 0.5 1.0 

Azimuth y [H] 

Figure 4. Flow pattern surrounding a protoplanet of dimensionless mass m = 10~ 2 in the linear regime, where perturbations are by 
construction small. Solid curves are contours of constant * = ^oo + ipyO, where tpyo is given by equation 1331 The far-field solution ('J'oo) 
does not include a headwind term (shear-only case where = 0). Pink circles indicate the perturbation in surface density: contours of 
<7 = 0.01 and 0.1 are shown. The Hill radius (r h = (m/3) 1 / 3 ) is given by the dashed circle and the Bondi radius (rt = m) is indicated 
by the black dot at the centre. These solutions are used as boundary condition for the numerical calculations in Section [4] 



in Section [T] is of special significance since it determines the 
horseshoe torque and connect to the planet's atmosphere. 
The horseshoe width x^ s is obtained by equating -i/>hs with 
the far-field solution, ^oo(xh s ), equation ([S]) : 



(35) 



(= 1.72ao\/q p /h in physical units; see Section [2.3[) . Equa- 
tion (|35p is in excellent agreeme nt with the findings of 
IPaardekooper fc Papaloizou |2009l ') who found a prefactor 



of 1.68. 

Figure [4] shows the flow pattern in the linearized ap- 
proximation with .Mhw = and m = 10~ 2 . Solid grey 
curves are iso-contours of SPoo + tp y o. The critical horse- 
shoe streamline $h s is highlighted blue. Iso-density lines, 
which follow from equation (|3ip . are in red. Figure 0] re- 
covers the global features of the perturbed shear flow (the 
horseshoe orbits in particular) at distances far from the 
planet (i.e., for r p£> m) as seen in many previous works 
(|Masset et al.ll2006l ; [Paardekooper fc Papaloizoull2009l ') . For 
large x the linearized approximation which relies on i < 
1 will break down. For even larger distances, i.e., for 

> 2/3, the PDE describing the system becomes hy- 
perbolic, signifying that a (density) wave is excited. Hy- 
perbolic (wave-li ke) equations require very d i fferent solu- 
tion techniques dGoldreich fc Tremainel Il98d : IWardl 1 19861 ; 
iTanaka et al]|200fl ). Including this feature is beyond the 
scope of this paper. 

For the purposes of this paper it is the flow structure 
in the vicinity of the perturber's Bondi radius that interests 
us. At these scales, the linear approximation is expected to 
break down and the nonlinear equations must be solved in 
earnest. The analytical solution in the linear regime pre- 
sented by equation (|33p will serve as the (outer) boundary 
conditions for the nonlinear calculation. 



4 STEADY-STATE FLOW SIMULATIONS 

In this section a strategy to numerically solve the general 
equations for S and * is discussed (Section 14. ip . Then, a 



parameter study is presented, varying the mass of the per- 
turber and the value of the headwind. 



4.1 Numerical algorithm 



Equations (1251) , (1261) are solved by iteration, considering one 
equation at a time. First an initial ^ is assumed, for ex- 
ample the linear solution, ^T/oo + tp y o- Equation ()26p is a 
scalar equation; given ^(x) and V\P it can be inverted to 
obtain E at every grid point. Then, using these values for E, 
equation (|25p - a diffusion equation - is solved for ^. This 
completes one cycle. 

Thus, when considering equation ([26} the RHS is 
known. The LHS depends on the model for the enthalpy W , 
for which an equation of state (EOS) needs to be specified. 
Obtaining E from equation ([26} is not without ambiguity, as 
is also illustrated in Fig. [21 Here the LHS of equation (126p . 
/ihs(E), is plotted as function of E for some arbitrarily- 
chosen value of IV^I 2 . Starting from E = 0, /m s (E) first 
decreases, reaches a minimum at (E m j n ,/ m i n ), and then in- 
creases towards infinity. Thus, for a given value of the RHS 
of equation (|26p there are either zero, one, or two solutions. 
In the case of two solutions, we always choose the largest E. 
This is simply because it gives the correct solution in the far 
field, provided no discontinuities (shocks) are present. 

However, during the iteration process, cases are encoun- 
tered in which the RHS of equation (|26p evaluates to less 
than /mi n . This is undesired. If this situation occurs, we re- 
adjust by 'mixing' it with the solution obtained in the 
previous iteration: 



* = w n <i> N + (1 - W„)*P, 



(36) 



with typ the previous solution for the stream function that 
did satisfy equation ([26}, ^jv the new solution, and w n a 
weight. Starting from w n = 1, w n is adjusted until ^ satisfies 
equation (|26p again. In runs where w n in this way is forced 
to become arbitrary low the iteration procedure clearly has 

failed. The calculations are then termi nated. 

The finite volume algorithm FiPy (|Guver et al.112009 1 ^ 



Downloadable at http://www.ctcms.nist.gov/fipy/ 
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Parameter 


oymbol 




Range 


Unit 


Resolution 


N ICS 


200 


100 — 400 




Headwind 






0, 0.05, 0.1 


Cs 


Mass 


m 


10" 2 


10" 3 — 10" 1 


h 3 M* 


Inner radius 


Ri 


0.1 


10" 3 —1 


Rb 


Domain size 


7" out 


0.3 


0.01—0.6 


H 


Transition mass 


St 


0.1 


10— 10 3 





Table 1. Default value and range of numerical and physical pa- 
rameters considered in this study. The last column denotes the 
unit to convert the dimcnsionlcss numbers into physical ones. 



is used to solve the partial differential equation (|25[) . A 
grid and boundary conditions (BCs) need to be specified. 
It is customary for these kind of simulations to consid- 
ered a rectangular grid. In this work, however, a polar grid 
(r = \J x 2 + y 2 \ 8 = arctan(y/a-)) is used with the perturb- 
ing body located at the centre of the coordinate system. The 
inner boundary of the grid is located at a radius n from the 
perturbing body and may correspond to the physical radius 
of the body. 

In FiPy either fr or the gradient 9 r fr must be specified 
as a BC. At the inner boundary (r = ri) the radial velocity, 
v r , vanishes. Therefore, dg^ — at the inner surface, and 
\& is constant (say fr = fri) on this surface. The value of 
fri is a priori unknown, but it can temporarily be set to 
an arbitrary value since equation (|25p will be unaffected by 
adding a constant to Therefore fri is fixed at *l/i =0. The 
outer boundary (r ou t) is then specified by Neumann BCs, 
i.e., dr^^out) = dri^ao + 4> v o) where the far field solution 
in polar coordinates reads froo = -|tOoor 2 cos 2 8 + v^r cos 6 
and ip y o is given by equation (|33ll . 

The obtained solution (with SET = 0) generally shows 
a mismatch at the outer boundary, i.e., 'l'(rout) ^ $ ou t, 
where the latter is given by the linear perturbation: 'tout = 
T ) + tpyo (r ou t ) . Let the difference between *]/ and frout 

be denoted C*. Although the constant CV is irrelevant for 
the velocity of the flow (which are derivatives of fr), it mat- 
ters for equation (|26|l , since it was assumed in its derivation 
that fr approaches fr<x> in the far field. The constant C* is 
obtained by minimizing the relative error at r — r out , 



E 



+ CV - fro 

*out.t 



(37) 



with respect to CV, where the index i refers to the i az- 
imuthal point on the outer surface. It is then added to the so- 
lution obtained by FiPy such that it approximately matches 
*out at r out . 



4.2 General results; key parameters 

By virtue of the various simplifications that have been em- 
ployed, the algorithm to solve for the steady-state flow pat- 
tern is computationally quite efficient; the CPU time is not 
much longer than a minute for a typical run on a standard 
desktop machine at the default resolution of 200 grid points 
in the radial and azimuthal direction. This allows us to carry 
out an extensive parameter study, studying the effects of 
the physical and parameter parameters see Table [1] The 
key physical parameters of interest here are: the headwind 



(fhw = A^hwCoo) and the planet mass (which is expressed in 
terms of the Hill radius ru, see Section ^, 31) . 

Not all runs converge. As expected from the discussion 
in the previous section, runs with a radial width that ap- 
proaches the scaleheight no longer converge, because the 
setup does not resolve shocks. Either one of equations (|25|) . 
(|26|) then fails to be satisfied. Typically, it is found that 
models with an outer radius (r out ) exceeding ~ 0.4 start to 
slow down significantly (because w n becomes low) and mod- 
els with r out > 0.6 fail to converge at all. The dependence 
on the outer radius is therefore investigated, see Section T4. 41 

The sensitivity of the flow pattern to the number of grid 
points in the radial and azimuthal directions (N r ,Ng) has 
also been investigated. Generally, we find that convergence 
is achieved by A r rca = 200, meaning that the N ICS = 400 
results did not provide a noticeably different result. But see 
Section f4. 41 for some exceptions. 

The flow and density structure interior to the Bondi ra- 
dius is affected by the inner radius Ri and transition density 
Et. These are also varied: Ri starting from the Bondi ra- 
dius until Rh/10 3 (which corresponds to the true physical 
radius at 5 AU); and Et is sampled at 10, 10 2 and 10 3 . The 
density close to R — Ri can become rather large; and due 
to the conservation of vortensity (a)/E), if non-zero, a very 
large azimuthal velocity emerges near the inner radius Ri. 
The large dynamic range in E and v in these cases often 
proves too much of a burden for the numerical code. How- 
ever, variations in Ri and Et are found not to significantly 
affect the flow pattern outside the Bondi radius. 



4.3 Runs without shear or rotation 
(headwind-only) 

In this section all parameters regarding rotation and shear 
are zero: cjoo = to — Q = 0. Consequently, the outer outer 
boundary condition fr ou t is set by the far field solution 
and does not involve ip y o, as the latter is only applicable for 
Keplerian shear. 



4-3.1 Incompressible runs 

When the flow is furthermore incompressible (E = E<x> = 1) 
equation (|25[) becomes Laplace's equation, V 2 fr = 0. For 
a cyl indrical geometry, the fl ow has an analytical solution 
(e.g., lLandau fc Lifshitjl 19591 ) : 



* = Mi 



2 i 

— ) cose = * c 

r 



A4hw^i cos# 



(38) 



where the unperturbed flow is assumed to move in the neg- 
ative y direction. This solution can be used to test the ac- 
curacy of the numerical model. 

Figure [5ji shows the stream pattern corresponding to 
the incompressible limit, where E is enforced to be unity. 
In the figure, the direction of the flow is indicated by black 
arrows, which indicate the direction of the flow and the flow 
velocity (larger arrows correspond to larger v but there is no 
absolute scale in any of the figures) . Contours of constant fr 
- streamlines - lie parallel to v. Note the vertical streamline 
at x — (9 = 7r/2) that must hit the objects due to sym- 
metry considerations. The surface potential (fri = fr(ri)) is 
therefore characterized by this value for fr; here, \&i = 0. 
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Figure 5. 2D inviscid flow calculations in the case of a pure headwind (no terms involving shear or rotation), (a) Incompressible flow 
with the Mach number of the headwind equal to 0.1; (b) gravitating flow with the Bondi radius equal to one-tenth of radius of the 
object; (c) gravitating flow at an increased Mach number of 0.5. In each panel the flow pattern is described by arrows and streamlines 
(gray curves) that are isocontours of <P. There is no absolute scale for the size of the velocity arrows among the panels. Solid, coloured 
contours indicate gas overdensities (E/Eoo) of: 1-2, 1.5 (salmon), 2.0 (thick red), 4.0 (salmon), 10, 10 2 , 10 3 , etc (all magenta). Note the 
focusing of streamlines starting in (b) and the distortion of the isodensity contours in (c). 



The computed flow pattern corresponds very well to the 
analytical solution (equation I38[) . but it does not give a per- 
fect match. The small offset in is caused by the boundary 
conditions at the (finite) outer radius. Here (at r = r ou t) 
d r ty = dr^oo, which deviates from equation (|38[) by a term 
of 0(r~ u 2 ). By increasing r ou t this error can be made arbi- 
trarily small, however. 

4-3.2 Compressible, gravitating flow 

In Fig. [Sp gravity has been including while keeping the flow 
irrotational. The value of the headwind («hw) now matters 
and it is fixed at 10% of the sound speed. The flow pattern 
is further determined by the ratio of the Bondi radius to 
the inner radius, Ri/Ri, = 0.1. Isocontours of E are nearly 
circular and displayed in Fig. [5p, c in red-shaded colours: 
E =1.2, 1.5, 4 (salmon), 2 (red), 10, 10 2 , 10 3 (magenta). 

For Rt, = R\ the flow does not significantly deviate from 
the incompressible limit. However, when the body increases 
in mass, it acquires a thick atmosphere. In the hydrostatic 
limit E(r) can be obtained by balancing the gravitational 
force to the pressure support: 



which results in E = Eoo exp(m/r) for an isothermal EOS. 
The radius r = m (which is the Bondi radius) signifies the 
point below which the gas density will increase exponen- 
tially. For r\ < m density gradients becomes very steep, 
and huge amounts of gas mass can be concentrated towards 
r = n as long as the EOS stays isothermal. Here however 
the EOS changes at a density Et = 10 2 (see equation I27p . 
causing a transition to a power-law dependence of E with r. 

As can been seen in Fig. [5p the (curl-free) flow focuses 
towards the gravitating body and slows down when it ap- 
proaches the body. This is simply a consequence of mass 
conservation: as E rises sharply, v has to decrease due to the 
steady condition that is imposed on the flow. As a result, the 



flow converges towards the centre, and then diverges again. 
The focusing of the gas gives rise to an hourglass pattern 
for the streamlines. 

For a Mach number of 0.1 in Fig. [5p, the density is 
well-approximated by the hydrostatic limit of solution (|39l) . 
Never the less the density isocontours are not entirely cir- 
cular, but the deviation is unnoticeable by eye. For larger 
values of TVfhw, however, the isocontours outside the Bondi 
sphere become more oval-shaped, see Fig. [SJ;, where the 
Mach number of the incident fl ow is increased to 0.5 . This 
trend agrees with the study of iLee fc StahleJ (|201ll ), who 
considered an 3D axisymmetric geometry. However, for cases 
where Mhw •C 1 the density can be well approximated as 
a function of radius only, E = E(r). This finding will be 
employed in Section [5] to find an analytical model for 



4.4 Shear-only runs; parameter study 

Next, the flow past bodies in a Keplerian potential is consid- 
ered: i.e., LJoo = —3/2 and uj = 1/2 (fio = 1) and the outer 
boundary condition at radius r out involves 4> y o (Section [3]). 
The headwind A'fhw = 0, such that the y-component of the 
flow velocity always vanishes at x = 0. The mass is taken to 
be m = 10 -2 , which implies a Hill radius of = 0.15. In 
physical units, this mass equals ~O.l-l.OM9 depending on 
the disc aspect ratio (see Fig. Q~} . The inner radius is again 
fixed at n = 0.1m (Ri = O.IRt), which implies an inner 
boundary somewhat larger than the physical radius of the 
protoplanet. The larger value is adopted for computational 
reasons; the value of n does not affect the flow pattern out- 
side the Bondi radius. 

In Fig. [6] the steady flow corresponding to these pa- 
rameters is presented at three different magnifications. In 
Fig. |6K the scale of the panel is that of the Hill radius, i.e., 
the range in both X and Y is 2Rh- For convenience, follow- 
ing Fig. [5] length units on both axis are given in Bondi radii. 
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X [R b ] X [R b ] X [R b ] 



Figure 6. Flow pattern and density contours in a shear-only run (neglecting the headwind contribution) for m = 10 -2 and an inner 
radius one-tenth of the Bondi radius. The panels give the flow pattern on scales of: (a) the Hill radius r^; (b) (c) . See 

Fig. [5] for a description of the contours. In addition, the solid blue curves correspond to ^ x'- the isocontour line of \P that intersects 
the stagnation point of the flow on the x-axis (where velocities vanish). This streamline defines the horseshoe region (where streamlines 
make U-turns) and the atmosphere of the planet (where streamlines are closed). 



Figure [Bp,c present zoom-ups of the flow in the vicinity of 
the planet. 

The topology of the flow features several qualitatively 
distinct regions. To the far-left and right, at scales \x\ > rh, 
it resembles the unperturbed solution: streamlines are nearly 
vertical with only little curvature. However, at smaller x- 
values the planet more strongly affects the flow, giving rise to 
the horseshoe region, quite similar to the linearized solution 
of Section [3] The streamline that divides these regions - the 
separatix streamline ^fx - is highlighted. At the point where 
it crosses the a-axis the flow stagnates: v = 0. The a;-value 
of this stagnation point is nonzero, in contrast to the linear 
solution. A new region, not present in the linear solution, 
therefore arises: the planet's atmosphere. 

Within the atmosphere region streamlines are closed 
and the gas, which encircles the planet, is b ound. The direc- 
tion o f the flow's rotation is prograde (e.g.. IrTAngero et al.l 
120021 ; iTanigawa et al.1 l2012h . For the shear only runs, it 
encompasses a region that is quite larger than the Bondi 
radius, r — m. Most of the atmosphere's volume - not 
necessarily its mass - is therefore of low (nebular) den- 
sity. At scales r < rb the density sharply increases. So 
does the velocity, which becomes more and more azimuthal 
(\v<t>/v r \ 3> 1). This increase follows from vortensity conser- 
vation, i.e., ui ~ wE. Nevertheless, the atmosphere at this 
stage is predominantly supported by pressure, instead of ro- 
tation. 

These res ults are br o adly in agreement with previous 
studies (e.g.. iMikil Il982l: iKorvcansky fc Papaloizoul Il99fj ; 



iBate et al.ll2003l ; iMachida et alj|2010h . although usually a 
more massive planet is considered (with a non-zero softening 
radius instead of a surface) and the focus lies on resolving the 
flow pattern on scales r S> H in order to resolve the spiral 
wave pattern excited by the Lindblad torques. Nevertheless, 
the general features de scribed above are all recovered. The 
steady state solution of IBate et al.l (|2003l ) assumes that the 
planet accretes mass in, essentially, the BHL-regime, which 
may be appropriate for large planet masses, but requires a 
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Figure 7. Dependence of the atmosphere size (the stagnation 
point where the critical streamlines intersect) i? a tm (given in 
terms of Bondi radius) on the numerical resolution N ICB (sym- 
bols) and the domain size r ou t (x-axis). Symbols that corre- 
sponds to the same r ou t are slightly offset for clarity. Results 



for the m 



10- 



runs lie above the dashed auxiliary line; runs 



with m = 10 — 2 below. This study indicates that convergence is 
achieved at a large domain size but that especially the low-m runs 
require a large resolution. 



very efficient cooling mechanism in the low planet regime 
(cf. discussion in Section [T]). 

The radius of the atmosphere JJatm may be defined by 
the point where the separatix streamlines intersect, i.e., it 
gives the distance to the stagnation point. Judging from 
Fig. [5] the atmosphere radius is ^3.5 Rb- In Fig. [7] the de- 
pendence of i? a tm against the two key numerical parameters 
- the numerical resolution N rcs and the domain size r out - 
is investigated both for an m — 10~ 2 and m = 10 -3 per- 
turber. The numerical resolution varies from N ICS — 100 
(these runs take a few seconds to complete on a desktop 
machine) to iV rcs = 800 (runs can approach one hour). The 
domain size r out goes up to r out = 0.6 scaleheights - the 
maximum where the runs still converge (see Section |4T}. In 
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Figure 8. Same as Fig.[6]but also including a headwind of i>i lw = O.lcoo {Mhw = 0.1). There are two solutions for the scparatix streamline 
(in blue): one associated with the horseshoe region, which lies at x ~ —10 in (a); and one circumscribing the planet's atmosphere. Note 
in (c) that the atmosphere is asymmetrical with nebula material penetrating within a fraction of the Bondi sphere. See the caption of 
Fig. [5] for the description and labeling of the contour levels. 



this figure, the atmosphere radius is also measured in Bondi 
radii. Note that the m = lCP' 5 runs (which lie above the 
dashed auxiliary line) have a Bondi radius that is a factor 
ten smaller than the m = 10~ 2 runs. Thus, in physical units 
R&tm is (as expected) larger for larger m. 

Figure [7] shows that in most - but not all - of the runs 
the numerical resolution hardly affects 7? a tm. More signif- 
icant is the dependence on the domain size r out . Figure [7] 
shows that a larger domain size results in larger atmosphere 
radii. Clearly, the outer boundary constraints influence the 
flow pattern to some extent and to be assured of conver- 
gence r out must be chosen as large as feasible. However, for 
the m = 10~ 3 runs the dynamic range, i.e., the ratio of r out 
to the Bondi radius (m) becomes very large in this way and 
a high resolution is needed to solve for the flow pattern at 
every point. But Fig. [JJ does show that all models converge 
for large r ou t and iV ros . In the following we will therefore 
stick with the standard of r ou t = 0.3 and iVr OS = 200 for the 
numerical parameters for which we expect the error in Tiatm 
to be at most 10%. 



4.5 Runs including shear and headwind 

Figure [8] shows the steady state flow pattern of a run that in- 
cludes a nonzero headwind (Mhw = 0.1) and a (Keplerian) 
shear term, again for a mass of m — 10~ 2 . The flow topol- 
ogy is quantitatively different from the shear-only case. The 
combination headwind+shear destroys the symmetry that 
was previously present along the x — axis; the stagnation 
points are no longer symmetrical around x = 0. As a result, 
there are two separatix streamlines and two solutions for 
one that encloses the planet's atmosphere ($ at m); and 
one, offset from the planet, that corresponds to the horse- 
shoe region (^hs). The planet's atmosphere is sandwiched 
by flow which moves in the direction of negative y - an ef- 
fect that is also seen for lower «t w . This topology is therefore 
the rule and the Uhw = run of Fig. [S] a special symmetric 
limit. 

The position of the separatix streamlines changes ac- 
cording to the planet mass and the value of the headwind. 



Figure [9] investigates the effects of varying the headwind 
(Alhw = 0, 0.05 and 0.1; from top to bottom) and the planet 
mass (m = 3 x 10" 3 , 10~ 2 and 3 x 10" 2 ; from left to right). 
Thus, Fig. [9^ corresponds to Fig. [6] and Fig. [9jd to Fig. [8] 
The scale of each panel is again the Hill radius. Note that 
the left-most highlighted streamline in the runs with a mod- 
est headwind (Mhvr = 0.05; panels d-f) is not associated 
with any stagnation point. Instead, it just happens to have 
the same value as that of the atmosphere streamline ^ a t m 
that encloses the atmosphere of the planet. 

Despite its importance, many works studying the flow 
behavior tend to ne glect the headwind term (but see, e.g., 
|Paardekoop"eill2009l for the implications of a large headwind 
on planet migration). In particular, the headwind may af- 
fect the accretion efficiency of small particles. Comparing the 
-Mhw 7^ to the Mhw = runs, it is noted that in the former 
the atmosphere region is much reduced. Small particles can 
be brought closer to the planet, even within a fraction of the 
Bondi radius. From this observation it seems that a (mod- 
est) headwind may actually promote protoplanet growth. 
Section [6] assesses this issue in more detail. 

The asymmetry of the atmosphere is quantified in 
Fig. 1101 which plots the radius of the atmosphere as func- 
tion of the mass of the planet for runs with and without a 
headwind. The size of the atmosphere R^ tln (R^tm) is de- 
fined to be the absolute value of the x-coordinate where 
f a tm intersects the positive (negative) z-axis. For A4hw = 0, 
-^atm = ^atm is seen to follow a trend that lies between the 
Bondi radius and the Hill radius. Thus, for a small planet 
it is significantly larger than the Bondi radius (cf. Fig. [6]). 
For Mh-w the symmetry breaking causes two solutions 
with i?^ m > R^tm- The negative solution always lies within 
the Bondi radius. In addition, R~ tnl and R^ tTn decrease with 
increasing headwind. 



5 ANALYTICAL DESCRIPTION OF THE 

FLOW PATTERN NEAR THE PERTURBER 

From the numerical results, it was found that the density 
E is primarily a function of r. This can be understood 
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Figure 9. Trends with increasing planet mass (from left to right) and increasing headwind (from top to bottom). Each panel covers a 
region two Hill radii in width, centred on the planet. 



from comparing the terms appearing in equation (|26[1 . When 
r < m the term representing the planet's potential, m/r, is 
dominant, whereas for large radius the first order solution 
of equation (|26[) is, by definition, "i/aa- At r ~ m, then, the 
Bernoulli equation approximately reduces to VK(E) = m/r, 
which in the isothermal regime leads to E = exp(m/r). 

The solution E = exp(m/r) corresponds to an isother- 
mal atmosphere in hydrostatic equilibrium. In the following, 
the hydrostatic solution is supposed to hold for the density 
structure, but not for the flow itself which is steady, but not 
static. The diffusion equation (|25|) is then left to solve. It 
can be written as 



-VlogE- V* + V 2 # = -E(&E - 2) (40) 



or in terms of the derivatives towards r and 9: 

r 2 6Vr*+ (r - r 2 dl ^ S ) dr^+dgg^ = -Er 2 (£E-2), (41) 

when the density is a function of radius only. The term in 
brackets on the LHS becomes (r+m)9 r \E' for E = exp(m/r). 

This partial differential equation (PDE) can be solved 
by adopting the Ansatz that is a superposition of a ho- 
mogeneous solution, $hom(r, 0), and a particular solution 
#r(r): 

*(r,0) = * hom (r,(9) + * r (r); (42) 
i.e., where $hom satisfies the homogeneous PDE 
r 2 d„-*hom + (m + r)a r *hom + <9 ee *hom = 0, (43) 
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Figure 10. Atmosphere size as defined by the distance to the 
stagnation points. The solution corresponding to the positive x- 
axis (-R^ tm ) is plotted by black symbols, while gray symbols cor- 
respond to itll m . In the shear-only case (A-I^w = 0) symmetry 
ensures that R^ tTn = K~ tm . A headwind breaks the symmetry 
and results in a noticeable decrease, of especially R~ tm - 



and ^r(r) the ordinary differential equation (ODE) 



r 2 *" + (m + r)^' r 



m/r 2 n 

e ' r 2 ■ 



~ m/r 
■ LOB ' 



(44) 



where primes denote derivatives to r. 

This simplifies the problem greatly. Equation (|43[) re- 
sembles Laplace equation but contains an additional d r ^ 
term. Equation (|44[) only contains r as the variable. Closed 
form solutions are possible. 



5.1 Solving equation (|43|) 

The strategy to solve this PDE is the same as with 
Laplace's equation. Assuming that the solution is separable, 
* hom (r,fl) = f(r)g(8), equation (|43)l reduces to two ODEs: 



r 2 /" + (m + r)f - n 2 f 
for the radial part; and 







d 2 g , 2 







(45) 



(46) 



for the angular part. Thus, g(9) — cosnO, where the sine so- 
lutions are discarded because of symmetry considerations - 
the far field solution ty^, only contains cosine terms. Clearly, 
n = ±1 correspond to the case of uniform flow (headwind) 
and n = ±2 holds for the shear case. For \n\ ^ 2 the solu- 
tions to equation ()45|) . denoted f„(r), read: 



For r»m, fn(r) is of 0([r/m] n ), also for negative n. The 
general solution for $hom therefore reads 



*hom(r,6i)= ^2 Cnf n {r)coan9, 



(48) 



with c n a constant that must be obtained from the boundary 
conditions. 



5.2 Solving equation (1441) 



Equation (|44[) is only of first order in the derivative of 
$[., and exhibits an analytical solution. For it reads 



K(r) = 



;^ [m 2 £Ei (f) - Cjr(m + r)e^ + 2(r 2 + Ci)j 



2r 



(49) 



where Ei(x) is the exponential integral (see equation IB2|) 
and C\ a constant of integration. This constant determines 
the direction of rotation (prograde or retrograde) of the gas 
flow within the atmosphere. Equation (|49|l can be integrated 
once more to provide (see Appendix [BJ, which introduces 
another constant of integration, Ci- 



5.3 Setting the integration constants 

We first consider the limit r 3> m, where the solution \f r hom + 
ty r should converge to the far-field solution, ^00- In the limit 
r>m, equation (|49[) can be approximated as 



3_ N 
2^ 1 ' 



independent of C\ . Integrating once more: 



(50) 



(51) 



where Wo, has been used instead of ui in the first term. Only 
terms of ©(r 1 ) or higher are included in the expansion. For 
r 3> m, the linear term is small and will be neglected. Sim- 
ilarly, the expansion of ^hom gives 



C2- 



■ cos 26> 



Ci — cos 9 

771 



(52) 



It is now required that the solution "thorn + matches the 
far- field solution "l/oo (equation [5]) , which can be written as 



■ COS 26 + «hwf cos ( 



(53) 



U(r) 



12 

2r 
— ( 

m 

1 



m/r 



2r 
m 



6r' 



V r I 1 



1 + — 



r 
in 
2r 
3m 



4r 

+ 1 + — 

m 



(>;- 



-2) 
-1) 



1 = 0) 
, = 1) 

i = 2) 
(47) 



The three terms correspond respectively to the n — 2 (shear) 
term of the homogeneous solution, the n = 1 (headwind) 
term, and the leading term of equation (|51[) . Thus the in- 
tegration constants of the homogeneous solution, ci and C2, 
simply follow from the requirement that the leading terms 
should match in the far field: 



C2 = 



(54) 
(55) 



The c_i and c_2 constants can be found by demanding 
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Figure 11. Flow patterns drawn from the analytical prescription for ^ outlined in Section [5] for perturber masses of: (a) m = 3x 10 — 3 , 
(b) m = 10~ 2 , and (c) m = 3 X 10 — 2 . The headwind is Ai = 0.05 in each panel and the control parameter y c is set at 0.05 as argued in 
the text. These panels should be compared to the middle row of Fig. [9] 



that Vl/ at the inner boundary (r = ri) is a constant, i.e., no 
azimuthal dependence. This gives 



C-2 = 



C-l = 



-C 2 



-Cl 



Mn) . 

/-a(n)' 
A(n) 



(56) 
(57) 



7-iW 

But these terms do not affect the solution for r ~ m and are 
rather unimportant. One may neglect these and the negative 
n terms of equation (]48p . 

Much more critical are the integration constants of 
Ci and C2. There is, unfortunately, no clear constraint that 
allows a straightforward determination. Somewhat more 
loosely, we can argue that the solution for ^ r + ^hom should 
match to the solution in the linear regime, 'I'oo + ip y o, at 
some point. In particular, let us consider the horseshoe re- 
gion for x — 0. Thus, we consider the y-axis (8 = n/2) and 
require that at a certain point (x, y) = (0, y c ) *I'hom+V'(j/c) = 



C 2 -^% +*r(2/c) =i>yo(yc 



(58) 



(Note that = for x = and cos27r/2 = — 1.) In addi- 
tion, it is required that the radial derivatives of the stream 
function match at this point: 



-2C2— \ + (9r*r)(y C ) = dyl]}yo(yc)- 



(59) 



Equation (|59p will determine Ci, whereas equation (|58|) 
determines Ci. The integration constants follow straight- 
forwardly from the ^ r and expressions, although the 
very expressions are rather cumbersome. Assuming that 
m < !/c C 1 simplifies these significantly, however, and 
we may approximate: 



Ci w my c 



log 



Vc \ 

V2J 



+ 7 



C 2 w Ci ( — + 7 ) + log (y c ) (m»c - Ci) 



(60) 
(61) 



im [(-2 + 2 7 - log(2))y c + >/2ir] + Ci log(m), 



where 7 is the Euler-Mascheroni constant. 

The choice for the integration constants has now been 
rewritten in terms of y c . Although there is no natural choice 
for y c , it is expected that the hydrostatic assumption for p 
will become less accurate for larger r. Judging from Fig. |3j 
a typical scale will be y ~ 0.1. Thus, y c may represent the 
point where the nonlinear solution transfers to the linear 
solution of Section [3] Nevertheless, there is some freedom 
allowed in choosing the precise value of y c . We find that 
with y c — 0.05 the solution very well matches those of the 
numerical simulations and we stick to this value. 

In summary, we have derived an analytical approxima- 
tion of the perturbed flow pattern for a low mass planet 
(m <C 1) under the assumptions that the density is ra- 
dial, consistent with the isothermal EOS, E = exp(m/r), 
and that its solution reaches that of the unperturbed flow 
at r> m. The perturbed stream function is equation (|42|l 
with equation (|51[l for ip(r) and equation (|48|l for ^hom with 
integration constants given by equations (|54|l - (|60p . 



5.4 Synthetic streamline patterns 

Figure [11] presents the streamline patterns derived from the 
analytical expressions presented above for the values of the 
mass parameter m. The free parameter y c was fixed at 0.05. 
This figure can be compared to the middle panels of Fig. [9] 
the corresponding numerical solutions. The analytic solution 
retrieves most features very well. The analytical description 
is advantageous, because it is much faster than running a 
full numerical simulation and does not require interpolation 
of the result. 



6 TRAJECTORIES OF SMALL PARTICLES 
AND THEIR ACCRETION BEHAVIOR 

Having quantified the velocity structure around the planet 
in the form of the analytic solution presented in Section [S] 
trajectories of solid particles can be calculated. Particle tra- 
jectories will deviate from those of the gas (the streamlines) 
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due to their inertia. That is, it takes a solid particle a char- 
acteristic time - the friction or stopping time t p - to align 
its motion with that of the gas. The dimensionless equiva- 
lent to t p is r p — tpQa- Big particles (t p 2> 1) like km-size 
planetesimals experience the gas only in the form of a small 
perturbation, which none the less is important to damp the 
random motions (eccentricity) of these bodies. Here, we con- 
sider small particles of r p -C 1. To first approximation, they 
tend to follow the gas. However, their not fully perfect cou- 
pling to the gas causes them to slowly settle due to the gravi- 
tatio nal force, e.g., t owards the central sta r (|Weidenschillina 
1 19771 ) or to planets (|Ormel fc Klahrll20ld ). 

The particle size r p adds another dimension to the pa- 
rameter space (along with m and A4hw) and it is left to 
a follow-up work to conduct a full exploration. Figure [12] 
only presents a preview, where particle trajectories - not 
gas streamlines - are plotted for m = 1CP 2 , .Mhw = and 
0.05, and r p = 10 -4 , 10 -3 , 10 -2 . These correspond to small 
particles in the ~100 ^m-10 mm size range, dependent on 
the values of the local gas density, disc radius, and particle 
internal densitjQ. 

4 The relation between the size s and friction time t p depends 
on the gas drag law. In the Epstein regime (small particle, low 
gas density), t p = pss/pCoo where s is the particle size. When 



In Fig. [T^] each curve represents a trajectory of a single 
particle. We have omitted drawing arrows to indicate their 
directions; but these are clear from comparison to Fig. [9p 
and e. Highlighted curves indicate particles that are accreted 
by the protoplanet. They all settle to the planet due to in- 
spiraling motion in the prograde manner. This means that 
the accretion rate only depends on the mass of the planet, 
not on its radius. This accretion mechanism qualitatively 
differs from the way planetesimals are accreted, where one 
relies on a sufficiently close encounter to 'hit' the target. 
Accretion through settling may provide a very fast channel 
for growth. For r v ~ 1 particles it has already been shown 
that their accretion radii are on the order of the Hill radius 
jOrmel fc KlahfeOirMLambrechts fc Johanserj|2012l ). How- 
ever, for smaller particles accretion cross sections are likely 
smaller since they couple more strongly to the gas. 

Indeed, Fig. 112b . d shows that the t p = 10~ 4 particles 
are difficult to accrete since their trajectories are strongly 
tied to the gas. The particle trajectories in these figures 
strongly resemble the gas streamlines of Fig. [9] In both 

particle sizes become larger than the mean-free-path i m f p of the 
gas, t p increases by a factor ~ Vmfp- For even larger particle 
sizes tp starts to depend on the particle-gas velocity itself; but 
this regime is not considered here. 
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Fig. 112b . d only two particle streams penetrate the planet's 
atmosphere and settle slowly towards the center (at some 
point when it is clear that the particle will be accreted, the 
integration is stopped). Most particles that approach the 
planet, however, cannot break through the atmosphere as 
they, being strongly coupled to the gas, are dragged along 
gas streamlines. 

IPerets fc Murrav-Clavl (|201ll ) introduced the wind- 
shearing radius i? ws to quantify the competing effects of 
the (central) gravitational two body force and the gas drag 
force (which pulls particle away from accretion). That is, 
only particles with impact parameters less than i? WB can be 
accreted. A similar wind-shearing radius can be defined for 
low-mass planets surrounded by an atmosphere. First, it is 
noted that only particles near the atmosphere boundary will 
be accreted. Therefore, particles traveling along the 'I'atm 
streamline are most likely to be accreted. The wind-shearing 
radius then becomes the bandwidth AJ? wa over which the 
gravitational pull of the planet is sufficient to drag them 
into the atmosphere region. This idea will be explored fur- 
ther in a follow-up work. 

Since the drag force equals F = (v p — v gas )/£ p , its influ- 
ence will decrease with increasing friction time. For particles 
of friction time t p = 1CP 3 (the middle panels of Fig. I12[) 
more particles accrete, but more so for the headwind case 
(7 counts) than the A4hw = case (4). Even for particle 
streams that do not accrete, deviations from the gas stream- 
lines can be noted in Fig. 112b and e. For the r p — 10~ 2 par- 
ticles in Fig. 112b and f the deviation from the streamlines 
is much more evident; close to the planet these particles de- 
couple from the gas and the impact parameter for accretion 
is very large. Note finally that due to the prograde rotation 
of the gas in the atmosphere, small particles that enter the 
atmosphere will also start to rotate in a prograde manner; 
and once they have settled towards the central object they 
will convey this angular momentum to the accreting body, 
causing it to spin up. Small-particle accretion therefore nat- 
urally explains the o bserved preference for prog rade spins of 
solar system bodies l|johansen fc Lacerdd 120101 ). 

Obviously, Fig.[l2]only presents a small subset of the pa- 
rameter space and the above argumentation is rather quali- 
tative. In a future work we will return to the question of the 
efficiency of small particle accretion and a provide a more 
quantitative analysis. 

7 DISCUSSION 

7.1 Caveats and extensions of the numerical 
method 

We have described a method to compute the gas flow 
around gravitating bodies in steady state. The idea of the 
adopted approach is to solve for the stream function 'J 
- a scalar quantity which fully describes the flow in 2D. 
The global flow pattern agrees well with previous numerical 
studies. However, the method employs several idealizations 
(the flow must be inviscid, 2D, barotropic, subsonic, etc), 
and it would be interesting to investigate whether these 
conditions can be relaxed. We remark that a supersonic 
expansion of the method has a lready been presented by 
iKorvcanskv fc Papaloizoul (|l996l ). which would allow us to 
explore particle accretion for more massive planets. 



Since the focus of this work lies on embedded plan- 
ets, an extension to three dimensions is perhaps the most 
obvious avenue. Indeed, 3D effects could have some ef- 
fects, e.g., with re gard to the width of the horseshoe region 
(|Bate et all 120031 '). Although the stream function formula- 
tion can be easily gen eralized to axisymmetric configurations 
(Lcc fc Stahlerll201ll ). it is not commonly used for truly 3D 
geometries, like the shear flow consider in this work. But if 
we pursue this path, the 3D equivalent to equation (0) reads: 

Ev = V* X Vg, (62) 

where g(x) is a second stream function. In 2D g = z re- 
covers equation ([7]). In 3D, constants of ^ and g consti- 
tute surfaces, and the intersections of two surfaces defines a 
streamline. Following the procedure similar to Section [2] we 
can obtain a set of diffusion equations, one for each compo- 
nent of the vorticity. Unfortunately, the diffusion equations 
are now functions of two scalar quantities (*? and g) and 
therefore require a much more complex procedure to solve, 
probably involving an iterative approach. 

7.2 Sensitivity to boundary conditions 

A quite general, and somewhat surprising, feature of the nu- 
merical model is its sensitivity to the outer boundary con- 
straints (BC). Here, we have used the linear solution pre- 
sented in Section [3] to provide the flow velocity at the outer 
boundary. This turns out rather well, since the flow pat- 
tern is not very sensitive to the choice of r ou t, although we 
found that the atmosphere region converges only for large 
''out (requiring high numerical resolution). However, if the 
flow quantities at the outer boundary is changed to the back- 
ground flow pattern, 'I'oo, omitting i/) H o _ which means there 
is that v x = by construction - it is found that the flow pat- 
tern is very different. In some cases even retrograde rotation 
around the perturber is observed. 

Since the linear solution that we presented in Fig. [4] 

agrees with previous wor ks (e.g., iMasset et al.l I2006I ; 

iPaardekooper fc Papaloizoul I2009I ). we have confidence in 
the robustness of the nonlinear calculations that we derived 
in the later sections. However, we have neglected possible 
contributions of th e spiral density wave, i.e., from the Lind- 
blad torque region. IPaardekooper fc Papaloizoul (|2009T ) con- 
sidered the effect of the Lindblad torque (which originates 
at distances x > H) on the flow in the vicinity of the proto- 
planet. They found that is was justified to neglect the Lind- 
blad torque in the limit of zero softening of the gravitational 
potential (like in this work). Nevertheless, it is worthwhile to 
include the effects from the line ar density wave theory (e.g., 
IWardlll997l ; ITanaka et all 120021 ) as a boundary constraint, 
and to assess how this affects the flow structure near the 
planet. 

7.3 Feasibility of the laminar approximation 

The laminar approximation - and the neglect of turbulence 
- is another caveat. One stability criterion is the Richardson 
number, which measures the stability of a stratified flow: 

R . n 2 _ 9(r)(dp\ /(dveV 
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where N 2 is the Brunt- Vaisala frequency. Richardson num- 
bers below a critical value, Ri cr it — 1 /4, indicate that the ki- 
netic energy term due to shear will outweigh the (stabilizing) 
term due to density stratification. The flow then overturns 
and becomes unstable. We have calculated the Richardson 
number for the simulations presented in Fig. [9] In the far 
field [r 3> 1) it was found that Ri < Ri cr it; but this is 
due to our initial setup, where Eoc is constant (allowing for 
gradients in Eoo and other background quantities [u>oo, Coa] 
would be another obvious extension of the method.) How- 
ever, Ri increases with decreasing r and closer to the planet 
Ri > Ricrit- 

Other instabilities may still operate though. For exam- 
ple, the flow within the atmosphere may be convectively 
unstable. Flow past bodies characterized by large Reynolds 
numbers will generally become turbulent in the wake of the 
body. Vortices may develop along the horseshoe streamline 
l|Koller et al.ll2003T ). Or the gas may just be simply unstable 
to start with, e.g., due t o the magneto-rotational instability 
(|Balbus fc Hawlevlll99lh . All these effects render the lami- 
nar approximation rather fragile. 

But it would be difficult, if not impossible, to present 
an analytic framework for non-steady flow. Our calculations 
provide the first order effect against which more specific sim- 
ulations can be contrasted. For example, Fig. [12] shows the 
first-order effect of the interaction of particle with gas. Tur- 
bulence will change the particle trajectories; but, although 
their individual trajectories can be very different, it is yet to 
be shown that turbulent motions will result in a significant 
deviation to the time-averaged accretion rate. 



7.4 Implications for protoplanet growth 

Notwithstanding the canonical scenario (|Safronov|[l969h . 
where protoplanetary cores are thought to be built up from 
big km-size (or larger) planetesimals, protoplanetary growth 
may in reality be driven by small particles, because: 

(i) planetesimals will grind themselves down, befo re an 
embryo can accrete them (|Kobavashi et al]|2010l . [201ll ) or be 
trapped in resonances (and su bsequently be ground down; 
IWeidenschilling fc Davisll 19851 ). 

(ii) (sub)mm-observations imply that a large reservoir of 
the dust is typically observed i n the ~ mm/cm-size range 
(e.g.. lAndrews fc Williamsll2005l ) . 

(iii) growth by (big) planetesimals is slow. This is due 
to the negative feedback of a growing embryo on the plan- 
etesimal population (i.e., it excites the m to large eccentric- 
ity, rendering grow th less efficient; e.g.. lKokubo fc Idalll998l ; 
iFortier et alTl2007l ). 

For these reasons, it is relevant to study in greater detail the 
interaction between small particle and embryo. In a follow- 
up work we will perform a quantitative (parameter) study 
to obtain the accretion potential of small particles. However, 
from Fig. [12] two trends become clear. First, very small par- 
ticles (i.e., dust) do not make formidable building blocks. 
They just follow the streamlines of the gas, and will only 
accrete on to the planet if the gas does not so as well. Sec- 
ondly, accret ion rates rise steeply wi t h par ticle size. In a 
recent study. lLambrechts fc Johansenl (|2012h showed that a 
putative core accretes r v ~ 0.1 pebbles very rapidly with 
impact radii approaching the Hill radius, consistent with 



the findings in Fig. I12fc and f (the scale of the panels is 
the Hill radius). That this pebble accretion scenario results 
in high accretion rates is now widely recognized and may 
serve as an alternative to (the much slower) planetesimal 
accretion. However, one must not forget that t p — 1 parti- 
cle also drift the fastest l|Weidenschilling| I1977T 1 ; they may 
be remov ed from the system be fore they encounter a pro- 
toplanet (|Kobavashi et al1l2010l ). In this regard, somewhat 
smaller particles may offe r a better channel for growth (see 
lOrmel fc Kobavashill2012l for an assessment of the combined 
effects that affect the accretion behavior of small particles) . 

Protoplanets can also grow by contracting their atmo- 
spheres, allowing more gas to become bound. In this study 
this effect was crudely reflected by the transition density 
parameter Et: the higher Et the further the isothermal 
regime penetrates, and the more massive the atmosphere. 
More realistically, the atmosphere mass Ma is a function of 
the (solid and gas) opacity, accretion rate, equation of state, 
etc, and is o btained by solving th e full equations of stel- 
lar structure |Mordasini et al.ll2012T ). These calculations are 
usually carried out in ID, assuming radial symmetry. But 
it is only through multidimensional (preferable 3D) studies 
that a true assessment of the boundary bet ween the atmo- 
spher e and nebula gas can be obtained (e.g.. iLissauer et al.l 
I2009T ). In this work we find that the headwind plays a crit- 
ical role. A headwind renders the flow - and therefore the 
atmosphere - asymmetrical, see e.g., Fig. [St, which causes 
nebular material to penetrate the atmosphere within a frac- 
tion of the Bondi radius. The consequences of the asymmetry 
and, in particular, the small atmosphere radius have yet to 
be addressed by protoplanet atmosphere models. We remark 
that, although the importance of the asymmetry increases 
with increasing Wh w , it is a general feature - seen also at low 
headwinds. 



8 SUMMARY AND CONCLUSIONS 

In this manuscript the inviscid, steady state solution of a 
subsonic flow past gravitating bodies was considered. First, 
in Section [31 an approximate but analytic solution was pre- 
sented. This solution only holds in the linear regime, that is, 
at distances much larger than the Bondi radius Rt- The lin- 
ear solution served as a boundary condition for the full (non- 
linear) numerical calculations (Section [3]). These showed 
the importance of the atmosphere region (a nonlinear phe- 
nomenon), where the flow curls around the planet in the 
prograde direction. Using the findings from the numerical 
simulations, we next constructed a more complete analyt- 
ical approximation of the flow pattern, which captures the 
dynamics of the atmosphere region (Section[5]). This solution 
very well matches the numerical findings and may be used 
in subsequent studies for which the flow pattern at scales of 
Rb is important. 

The adopted framework is quite general, with param- 
eters describing the mass of the perturber, the headwind 
velocity «hw, and the shear parameter Uoo. Apart from Sec- 
tion 14.31 this study has focused on the flow pattern past 
small planets embedded in Keplerian discs = — 3Slo/2). 
However, flow patterns around 'bodies' occur in diverse as- 
trophysical settings, e.g., in (small) satellites e mbedded cir- 
cumplanetary disks (e.g.. [Estrada et al.ll2009l '). star forma- 



© 0000 RAS, MNRAS 000, 000-000 



The steady-state flow pattern past gravitating bodies 17 



tion jKrumhoIz et al.ll2005l). o r accretion disks around black 
holes ( McKernan et al.l 1201 ih . As long as these embedded 
objects do not open gaps and the system is in steady-state, 
the framework constructed in this paper also provides a de- 
scription for the flow pattern in these environments. 
The key findings of this study are the following: 

(i) For flows without shear, the steady-state flow pattern 
resembles that of an hourglass (Fig. [5]), where gravity focuses 
the streamlines. 

(ii) For Keplerian shear flows the vorticity, which is am- 
plified near the planet due to conservation of vortensity, 
plays the dominant role. There is a point in the simulation 
where the flow stagnates. In general there are two solutions 
for the associated stagnation streamline - one defining the 
atmosphere region and one defining the horseshoe region - 
which only coincide when «h w = 0. 

(iii) For the (general) Vh-w 7^ case, the boundary of the 
atmosphere is asymmetric and lies within the Bondi sphere 
(Fig. llQ|) . The size of the atmosphere reduces with increasing 
headwind but its averaged radius is best approximated by 
the Bondi radius Rt, (rather than the Hill radius). 

(iv) When the flow is subsonic (in particularly, w^w •C 1), 
the density near the Bondi radius can be well approximated 
by the hydrostatic solution, such that E becomes a function 
of radius only. In this approximation, an analytic solution 
for the stream function has been obtained (see Section [5] 
and Appendix [Bj . This solution can be used to compute 
the gas drag force which solid particles experience during 
their encounter with the planet. 

(v) When they can penetrate the atmosphere region, par- 
ticles can be accreted as they settle to the planet through 
the circumplanetary disc. Accretion of small particles (dust) 
is suppressed as they strongly couple to the gas flow, but 
we find that for ~mm-size and larger particles accretion 
rates should be high. There is also tentative evidence that 
a (moderate) headwind accelerates particle sweepup, due to 
the smaller atmosphere. 
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APPENDIX A: DERIVATION OF EQUATIONS 
(29), (30) 



The linearization of equation ([25} proceeds as follows. The 
term V*/E is replaced by V(*oo + V)/(l + °0 ~ W" + 
(1 — <r)V5 , 00 , where it is assumed that cr < 1 and that 
we can neglect terms like aVip. Since V^oo = (—oJooX + 
Mhw)ex only has an x-component, the divergence operator 
V- is replaced by a differential d x . In this way equation Q29p 
becomes 

X7 2 ip — (d x cr)(uJooX — Mhw) — ojoo(l — a) = —uioo — Cjcf. (Al) 

from which, after re-arranging, equation (|29p is retrieved. 

For the linearization of equation Q26| l. an isothermal 
EOS is assumed, such that W^E) = logE ss a. The kinetic 
term |V^| 2 /2£ is linearized as 

(A2) 

On the RHS of equation ((26} * = vE^ After noting that 
the first order terms cancel (they obey the far-field solution), 
equation (|30[) is obtained. 



APPENDIX B: FULL, NONLINEAR 
EXPRESSIONS FOR THE FLOW PATTERN 

In Section [5] an approximate solution for the flow in terms 
of the stream function (see equation [7]) was presented in 
terms of an homogeneous solution, *hom, and a particular 



solution, *(r): *(r, 6) = * h om(r, 0) + *(r). We found that 
(equation I49[): 



*'(r) = 



[m 2 £Ei (f ) - Cbr{m + r)e* + 2(r 2 + Ci)] 



2r 



(Bl) 



where C\ is an integration constant, and Ei(x) the exponen- 
tial integral 



Ei(x) = / dt — . 

J — CXI ' 

Equation (|B1|I can be integrated once more 
* r ( r ) = I [-Ei (™) (4Ci + m 2 £Ei (-) + 2m 2 

+4C 2 + 8m 2 £Ei ' 



(B2) 



+re m/r (2(m + r) - we m/r (4m + r))] 



(B3) 



These expressions contain two integration constant, Ci 
and Ci. As argued in Section these are chosen such that 
they match the linear solution at a point (x,y) = (0, y c ) or 
r = y c and 6 = -n/2 in radial coordinates. This resulted in 
conditions given by equations (|58|l . (|59|l . The derivative of 
equation (|33p . ip y o reads 



dytpyo — rn 



cos y Ci(y) - - siny (w - 2 Si(y)) 



(B4) 



where y = y/y/2, and y > and a Keplerian disk (uj = 1/2) 
have been assumed. 

The constants C\,C2 follow by rearranging equations 
(|58[1 . (|59|l . i.e., by isolating Ci and C2. This is a straight- 
forward procedure, but rather cumbersome. For a Keplerian 
disk (uj = 1/2) we obtain in this way 



C\ = e yc I my c Ci (y c ) cos y c — e vc 



my 



c [| - Si (y c 



1 2„. ( m \ 2 



sin y c 



1 2m 1 ] 

+^y c e yc (y c + m) + -y c (3y c + m) > . (B5) 

An expression for C2 can be obtained in a similar way. In the 
limit of m < |/ c < 1 these expressions simplify significantly 
and read: 



Ci w myc (log U^j + 7 
C* 2 « -^m 1 12\/27r + 4y c 



log I — ) + 67 - 6 



(B6) 
(B7) 



+ lo; 



+ 67 (- log (y c ) + log(m) + 7) 



8 

where 7 ~ 0.577 is the Euler-Mascheroni constant. 

To conclude, let us summarize the expressions for the 
flow velocity that follow from the full nonlinear model. For 
the radial component it reads: 

1 1 x ,t, ma\ 

V r = — = — 3fl*hom, (B8) 

2j r Lr 

where S = exp(m/r), and ^bom is given by equation ((48} 
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which includes the ci,C2 integration constants. The az- 
imuthal flow velocity reads 

V 9 = -^d r V = -~ [9r*hom + *»] (B9) 

and includes in addition the integration constant Ci through 
the fy'r term (equation IB 

This paper has been typeset from a TfrjX/ KTJ5X file prepared 
by the author. 
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